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The  distributions  of  potential  and  reaction  rates  in  a  lithium  ion  battery  during  discharge  process  have 
great  influences  on  the  battery  thermal  characteristics.  A  two-dimensional  electrochemical— thermal 
model  has  been  developed  for  a  cylindrical  LiFeP04  battery  by  coupling  the  mass,  charge,  and  energy 
conservations  as  well  as  the  cell  electrochemical  kinetics.  The  model  also  includes  battery  current  col¬ 
lecting  tabs.  The  modeling  results  are  validated  for  both  the  electrochemical  performance  and  thermal 
behavior  during  galvanostatic  discharge  process.  The  modeling  results  agree  well  with  the  experimental 
data.  The  placement  of  the  positive  and  the  negative  current  collecting  tabs  on  the  cylindrical  battery  was 
found  to  have  a  significant  effect  on  the  distributions  of  its  potential  and  local  reaction  rates,  which 
therefore  affect  heat  generation  rates,  and  thus  the  distribution  of  the  temperature  within  the  battery. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Large  capacity  or  large  size  lithium  ion  power  batteries  are 
widely  used  in  the  hybrid  electrical  vehicles  (HEV)  and  electrical 
vehicles  (EV).  The  performance  of  a  high  capacity  lithium  ion  power 
battery  is  influenced  by  the  geometrical  characteristics,  the  placing 
of  the  current  collecting  tabs,  and  the  local  reaction  rates  [1  ].  Larger 
geometry  sizes,  higher  heat  generation  rate  and  extreme  operating 
conditions  could  lead  to  quite  non-uniform  thermal  and  electrical 
behaviors  [1,2].  Therefore  large  capacity  lithium  ion  power 
batteries  require  a  more  comprehensive  understanding  of  their 
electrochemical  fundamentals  and  the  thermal  behavior. 

Modeling  and  simulation  of  lithium  power  batteries  is  becoming 
increasingly  important  in  developing  the  fundamental  under¬ 
standing  of  electrochemical  performances  and  thermal 
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characteristic  for  improving  battery  thermal  management  systems. 
In  recent  published  works,  several  one  dimensional  (ID)  electro¬ 
chemical-thermal  coupled  models  for  lithium  ion  batteries  have 
been  used.  These  models  comprising  conservation  of  mass,  charge, 
energy,  and  electrochemical  reaction  are  developed  based  on  the 
kinetics  models  [3-12  .They  usually  couple  the  ID  kinetic  model  of 
local  cell  unit  to  the  2D  or  3D  thermal  model  of  the  whole  battery 
cell  by  adding  heat  source  terms  and  temperature  dependent  pa¬ 
rameters,  without  considering  the  current  collecting  tabs  [7-10]. 
However,  little  attention  is  paid  to  the  modeling  and  simulating  of 
the  whole  battery  including  both  the  current  collecting  tabs  and  the 
sandwich  structures  of  cell  units,  which  are  closer  to  the  real  bat¬ 
tery  and  more  suitable  for  the  large  capacity  or  large  size  lithium 
ion  power  batteries  [12-15].  In  addition,  for  batteries  with  LiFePOq 
electrodes,  only  a  handful  of  kinetic  modeling  studies  have  been 
published  so  far  [8,12,11,16]. 

In  this  work,  a  two-dimensional  electrochemical-thermal 
model  has  been  developed  for  a  commercial  Type  38120  cylindrical 
LiFePCH  power  battery  by  coupling  mass,  charge,  and  energy 
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conservations  as  well  as  electrochemical  kinetics.  The  current  col¬ 
lecting  tabs  of  the  cylindrical  battery  are  considered  in  this  model. 
Effects  of  the  current  collecting  tabs  placing  on  the  cell  perfor¬ 
mance  are  also  investigated. 


2.  Coupled  electrochemical-thermal  model 

2.1.  Computational  domain  and  model  assumptions 

The  two-dimensional  electrochemical-thermal  model  for  a 
commercial  Type  38120  cylindrical  LiFeP04/Carbon  power  battery 
is  developed  based  on  the  principles  of  the  mass,  charge,  and  en¬ 
ergy  conservations  as  well  as  the  electrochemical  kinetics.  Table  1 
lists  the  application  of  these  conservation  equations  in  each 
respective  domain.  The  details  of  these  equations  will  be  given  in 
the  later  sections. 

The  nominal  voltage  and  capacity  of  Type  38120  are  3.2  V  and 
10  Ah.  Fig.  1  illustrates  the  two-dimensional  computational  domain 
of  the  cylindrical  battery  in  the  axial  cross-sectional  direction.  The 
domain  includes  several  different  regions:  the  cell  core  part  (battery 
jelly-roll)  with  the  structure  of  spiral-wound  cell  unit;  the  positive 
and  negative  terminals  as  well  as  the  tab  groups;  the  battery  outer 
can  made  of  stainless  steel.  The  spirally  coiled  cell  unit  is  wound  up 
in  the  form  of  a  jelly-roll,  so  it  appears  as  layered  structure  in  the 
two-dimensional  rectangular  axial-section  as  shown  in  Fig.  1(a),  (b). 

As  shown  in  Fig.  1,  the  basic  cell  unit  is  assumed  to  be  a  sandwich 
structure  that  is  composed  of  the  negative  electrode  (ne),  the 
separator  (sp),  the  positive  electrode  (pe),  and  the  two  current 
collectors  (cc)  including  the  tabs.  Each  positive/negative  tab  is 
welded  to  the  corresponding  current  collector  and  terminal  of  each 
basic  cell  unit.  The  electrodes  are  porous  solid  matrix  that  consists 
of  active  particles  with  spherical  shapes  of  uniform  sizes  and  ad¬ 
ditives  (Fig.  1(c)).  The  positive  electrode  contains  active  material 
particles  of  iron  phosphate  (LiFeP04)  and  the  negative  electrode 
contains  the  active  material  particles  of  amorphous  carbon  (LiCe). 
The  separator  is  a  porous  polymer  membrane  that  constitutes  a 
physical  barrier  between  the  two  electrodes.  Both  electrodes  and 
separator  porous  matrix  are  impregnated  with  electrolyte  ensuring 
the  transfer  of  lithium  ions  between  two  electrodes.  The  electrolyte 
is  a  mixture  of  cyclic  and  linear  carbon  solvents  in  which  a  lithium 
salt,  lithium  hexafluorophosphate  LiPF6  (1.5  mol  L-1),  is  dissolved 
in  a  nonaqueous  1:2  liquid  mixture  of  ethylene  carbonate  (EC)  and 
dimethyl  carbonate  (DMC).  The  electrochemical  reactions  that 
occur  at  the  interface  of  the  electrode  and  electrolyte  both  in  the 
negative  and  positive  during  discharge/charge  are  then, 

discharge 

Negative  Electrode  :  LixC6  Lix_zC6  +  zLi+  +  ze~ 

charge 

charge 

Positive  Electrode  :  LivFeP04  *±  Liv_zFeP04  +  zLi+  +ze 

discharge 

where  x  is  the  stoichiometric  coefficient  or  the  number  of  moles  of 
lithium  present  in  the  graphite  structure  (Ce),y  is  the  stoichiometric 
coefficient  or  the  number  of  moles  of  lithium  in  the  olivine  structure 


Table  1 

Governing  equations  of  different  subdomains. 


Governing  equations 

cc 

pe/ne 

sp 

Terminal/tab 

Outer  can 

Electrochemical  kinetics 

- 

V 

- 

- 

- 

Charge  conservation 

V 

V 

V 

V 

V 

Mass  conservation 

- 

V 

V 

- 

- 

Energy  balance 

V 

V 

V 

V 

V 

of  iron  phosphate  (FeP04),  Li+  is  the  lithium  ion,  z  is  the  number  of 
moles  of  lithium  taking  part  in  the  electrochemical  reaction. 

2.2.  Electrochemical  kinetics  at  the  interface 

The  local  charge  transfer  current  density  is  given  by  the  Butler- 
Volmer  equation  for  electrode  kinetics 

Jn=jo{exp(^??)-exp(-^??)}  (1) 

where  jo  is  the  exchange  current  density,  aa  and  ac  are  the  anodic 
and  cathodic  charge  transfer  coefficients,  rj  is  the  local  surface 
overpotential  and  F  is  the  Faraday  constant.  The  exchange  current 
density  is  given  by 

Jo  =  F/<0c“s  (c1>max  -  cljSurf)“ac“'surf  (2) 

where  k o  is  the  reaction  rate  constant,  Ci>max  is  the  maximum 
lithium  concentration  in  the  active  electrodes  and  c\  is  the  lithium 
concentration  on  the  surface  of  the  active  particles.  The  subscripts  1 
and  2  represents  the  solid  phase  and  electrolyte  phase  respectively. 
The  overpotential  is  defined  as 

V  =  4>\  —  4>2  ~  ^eq  (3) 

where  Ueq  is  the  open  circuit  potential  of  the  electrode,  0i  is  the 
solid  phase  potential  and  02  is  the  liquid  phase  potential. 

2.3.  Charge  conservation 


The  governing  equation  for  charge  conservation  in  the  positive/ 
negative  electrode  is  expressed  as 


V  •  ij  +  V  •  1*2  =  0 

(4) 

V-ij  =  -Sajn 

(5) 

V-l2  =  Sajn 

(6) 

where  tj  refers  to  the  electronic  current  density  in  the  solid  phase, 
i2  is  the  ionic  current  density  in  the  liquid  phase,  and  Sa  is  the 
specific  surface  area. 

2.3.1.  Electron  transport  in  solid  phase 

The  transport  of  electrons  in  the  solid  phase  follows  Ohm’s  law 
which  can  be  expressed  as  follows 

h  =  i  (7) 

where  refers  to  the  effective  electrical  conductivity  of  solid 
phase. 

2.3.2.  Lithium  ion  transport  in  solution  phase 

The  transport  of  lithium  ions  in  electrolyte  can  be  expressed  as 
follows  [3,17]: 

h  =  -afv<t>2  +  2Rf2  (l  +  0  -  4)V(ln  c2)  (8) 

where  h  consists  of  two  terms,  the  first  term  following  Ohm’s  law 
and  the  second  accounting  for  the  ionic  concentration  gradient,  j± 
is  the  average  molar  activity  coefficient,  t+  is  the  transferring 
number  of  lithium  ions  in  the  liquid  phase,  and  C2  is  the  lithium  ion 
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Fig.  1.  Schematic  graph  of  (a)  axial  cross-section  of  Type  38120  battery  (b)  partial  enlarged  detail  of  two-dimensional  computational  domain  showing  various  functional  layers  in 
the  jelly  roll  (c)  diffusion  of  lithium  (Li+)  in  active  material  particles  and  electrolyte  in  the  electrodes. 


concentration  in  the  liquid  phase,  R  is  the  universal  gas  constant, 
and  T  is  the  temperature  of  electrolyte. 


2.4.  Mass  conservation 


2.4 A.  Lithium  in  active  particles 

The  mass  conservation  of  lithium  in  the  spherical  active  material 
particle  can  be  described  by  Fields  law 


9Cj  1  d 
dt  +  r2  dr 


-rzDi 


=  0 


(9) 


where  C\  is  the  concentration  of  lithium  in  the  active  material 
particles  of  the  electrode,  t  is  the  time,  Di  represents  the  diffusion 
coefficient  of  lithium  in  the  solid  phase,  and  r  is  the  radial  coordi¬ 
nate  inside  a  spherical  particle,  which  is  no  larger  than  the  particle 
radius  R. 


2.4.2.  Lithium  ions  in  the  electrolyte 

The  mass  conservation  of  lithium  ions  in  the  electrolyte  can  be 
stated  as 


dc2  ^  j  5ajn 

e2-^  +  V-J2=-r  (10) 

with 

h  =  -Dfvc2+^t  (11) 

where  J2  is  the  molar  flux  of  lithium  ions  that  consists  of  two  terms: 
the  first  term  following  Fields  law  and  the  second  accounting  for  the 
electromigration,  represents  the  effective  diffusion  coefficient 
of  lithium  ions  in  the  electrolyte,  and  £2  is  the  volume  fraction  of  the 
liquid  phase.  The  term  on  the  right  side  of  the  Eq.  (10)  is  zero  when 
applied  in  the  separator. 


2.5.  Energy  balance 

The  energy  balance  in  the  lithium  ion  battery  can  be  expressed 
as  follows  in  which  there  are  three  parts  of  heat  generation  sources: 
reaction  heat  qrea,  activation  polarization  heat  qact,  and  ohmic  heat 

Qohm- 


a  t  7 

Pc P at  —  ^  T  =  Qrea  +  Qact  +  Qohm  (12) 

The  heat  generation  due  to  electrochemical  reaction  is 

<Jrea  =  =  Sa/nT^  (13) 

The  active  polarization  heat  is 

Qact  —  SajnV  (14) 

The  heat  generation  due  to  ohmic  losses  consists  of  electrical 
ohmic  heat  in  the  solid  phase  and  ionic  ohmic  heat  in  the  liquid 
phase 

Qohm  =  -ir'Wq  -12-V02  (15) 


2.6.  Temperature  and  concentration  dependent  parameters 

The  positive  and  negative  electrode  active  materials  are  Liy_ 
FeP04  and  Li*C6  respectively.  The  electrolyte  is  1.5  mol  L-1  LiPF6_ 
dissolved  in  a  mixture  of  2:1  EC/DMC.  The  model  input  parameters 
are  listed  in  Table  2  including  battery  geometrical  and  design  pa¬ 
rameters,  kinetic  and  transport  properties  and  thermal  properties. 

2.6 A.  Electrode  kinetics  related  parameters 

2.6  A  A.  Reaction  rate  constant  k0.  The  temperature  dependence  of 
reaction  rate  constant  k0  follows  the  Arrhenius  formula  [4,8,9,13]: 
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Table  2 

Model  parameters  for  a  10  Ah  LiFeP04  power  battery. 


Parameter 

Unit 

Al  foil 

Cathode 

Separator 

Anode 

Cu  foil 

Outer  can 

Ref. 

e\ 

- 

0.435 

- 

0.56 

- 

- 

Estimation 

e 2 

- 

0.306 

0.54 

0.268 

- 

- 

Estimation 

<5; 

|im 

20 

91 

40 

142 

20 

Estimation 

r0 

|im 

- 

1.15 

- 

14.75 

- 

- 

[11] 

t-2,0 

mol  m  3 

- 

1500 

- 

- 

[19] 

Cl.o 

mol  m  3 

- 

3900 

- 

16,361 

- 

- 

[19] 

Cl, max 

mol  m  3 

- 

22,806 

- 

31,370 

- 

- 

[19] 

aa 

- 

0.5 

- 

0.5 

- 

- 

[11] 

ac 

- 

0.5 

- 

0.5 

- 

- 

[11] 

D ! 

m2  s"1 

- 

1.18  x  10  18 

- 

3.9  x  lO"14 

- 

- 

[9,19] 

D2 

m2  s"1 

- 

Eq.  (21) 

- 

- 

[10] 

EaR 

kj  mol-1 

- 

4 

- 

4 

- 

- 

[11] 

EaD 

kj  mol  1 

- 

20 

- 

4 

- 

- 

[11] 

(T\ 

S  m1 

3.83  x  107 

0.1 

- 

100 

6.33  x  107 

- 

[1] 

02 

S  m1 

- 

Eq-  (19) 

- 

- 

[1] 

ko 

m2-5  mol-0  5  s"1 

- 

3.626  x  10"11 

- 

1.764  x  10"11 

- 

- 

[9] 

4 

0.363 

[9] 

k 

W  m  1  K1 

1 60 

1.48 

1 

1.04 

400 

44.5 

[9] 

P 

kg  m"3 

2700 

1500 

492 

2660 

8700 

7850 

[11] 

4,i 

J  kg"1  K-1 

900 

1260 

700 

1437 

385 

475 

[11] 

CP,2 

J  kg"1  K-1 

- 

2055 

- 

- 

[11] 

Tref 

K 

298.15 

F 

C  mol-1 

- 

96,487 

- 

- 

[11] 

I<0(T)  =  feo.ref  exp 


EaR 

R 


(16) 


where  EaR  is  the  reaction  activation  energy,  and  /<o,ref  is  the  reaction 
rate  constant  at  reference  temperature  Tre f. 


2.6A.2.  Open  circuit  potential  (OCP)  Ueq.  The  open  circuit  potential 
(OCP)  depends  on  both  temperature  and  concentration,  which  can 
be  expressed  with  a  Taylor’s  series  expansion  [13]: 

Ueq  =  Ueq,ref+^(r-rref) 

=  ueq  ref  (soc,  Tref)  +  (r  -  Tref )  (17) 

where  the  state  of  charge  (SOC)  of  both  positive  and  negative 
electrodes  can  be  determined  as  follows: 

SOC  =  SL-  (18) 

C1  ,max 

As  shown  in  Eq.  (17),  both  the  open  circuit  potential  (OCP)  and 
the  entropy  change  AS  depend  on  the  state  of  charge  (SOC).  The 
dependence  of  these  two  variables  on  the  SOC  could  be  determined 
by  experiments.  Fig.  2  represents  AS  and  OCP  changing  curves  with 
the  SOC  at  the  reference  temperature  25  °C  for  LiyFeP04  and  LixC6 
electrodes,  respectively. 


(,+ffi|)(,-,*)-F!ili(0-60,-a24>‘00,‘?!+0-982 

x  [1-0.0052(7-294)] 
x  (0.001c2)15) 

(20) 

The  Li+  diffusion  coefficient  in  the  liquid  phase  D2  is  given  in  Eq. 
(21)  [10,20-22]: 

D2(C2J)  =  1 0  (~4-43-  (’  So  001C2)  -0-22  X0.001C2)  -4  ^ 


2.6.3.  Solid  phase  parameters 

2.6.3 A.  Li  diffusion  coefficient  Dj.  The  temperature  dependence  of 
the  Li+  diffusion  coefficient  Di  in  the  solid  phase  follows  the 
Arrhenius  formula  [4,8,9,13].  The  constant  Di>ref  at  the  reference 
temperature  Tref  is  listed  in  Table  2. 


Di(T) 


Dl.ref  exp 


Eap 

R 


(22) 


where  FaD  is  the  activation  energy  for  diffusion. 


2.7.  Boundary  and  initial  conditions 


2.6.2.  Liquid  phase  parameters 

Eqs.  (19)  and  (20)  describe  the  temperature  and  concentration 
dependence  of  ionic  conductivity^  and  thermodynamic  factor 
(1  +  aln/±/alnc2)(l  -  t+)  for  LiPF6  in  EC/DMC  (2:1  by  volume) 

[10,20,21]: 

<72(c2,7)  =10 -4c2  x  1.2544  x  (-8.2488  +  0.0532487 

— 0.0000298772  +  0.26235  x  0.001c2  -  0.0093063 
x0.001c27  +  0.000008069  x  0.001c272  +  0.22002 

xl0~6c2  —0.0001765  x  l(T6d;7)2 

(19) 


At  the  electrode/separator  interface  I,  insulation  is  specified  for 
the  electronic  current  of  the  solid  phase.  Continuity  is  specified  for 
the  ionic  current,  lithium  ion  flux  of  liquid  phase  and  heat  flux. 

I :  n  h  =  0,  rw'2|I+  =  n-i2|,_,  n-J2|I+  =  n-J2 nq|I+  =  n-q|,_ 

(23) 

At  the  electrode/current  collector  interface  II,  insulation  is 
specified  for  the  ionic  current  and  lithium  ion  flux  of  the  liquid 
phase.  Continuity  is  specified  for  the  electronic  current  and  heat 
flux. 

II  :n  *2  =  0,  n-!'1|n+  =  n-i1|[I_,  n-J2  =  0,  n  q|n+  =  n-q|„_  (24) 

At  the  battery  jelly-roll/can  interface  V,  continuity  is  specified 
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Fig.  2.  Thermodynamic  properties  used  in  the  model  as  a  function  of  SOC  in  the 
electrodes:  (a)  entropy  change  [5,18],  (b)  reference  open  circuit  potential  [19]. 


for  the  heat  flux.  Insulation  is  specified  for  the  lithium  ion  flux, 
electronic  current  and  ionic  current. 

V :  ni2  =  n  il  =  n-J2  =  0,  n  q|v+  =  nq|v_  (25) 

At  the  current  collector/tab/terminal  interface  VI,  insulation  is 
specified  for  the  ionic  current  and  lithium  ion  flux  of  liquid  phase. 
Continuity  is  specified  for  the  electronic  current  and  heat  flux. 

VI :  n  •  i2  =  n-J2  =  0,  n  i!  |VI+  =  n-i,  |V|_,  n  <?lvi+  =  «  Qlvi 

(26) 

At  the  positive  terminal  boundary  III,  the  applied  current  of  the 
battery  is  defined  as  follows 

III  :  n  It  =  I'app  (27) 

The  negative  terminal  boundary  IV  is  grounded,  thus 

IV  :  0!|IV  *  0  (28) 

The  working  voltage  of  the  battery  can  be  defined  as 

^cell  =  01  IlII  -  01 1 IV  (29) 


V:-A-  =  h(T-Tamb)  (30) 

where  the  lumped  heat  transfer  coefficient  h  is  determined  by 
comparing  experimental  data  of  the  battery  average  surface  tem¬ 
perature  with  simulated  results  in  Fig.  3.  A  lithium  ion  battery  was 
placed  in  a  thermostatic  chamber  at  55°  Celsius  for  10  h,  ensuring 
that  the  overall  temperature  of  the  battery  is  uniform.  Then  the 
battery  was  put  in  the  environment  with  the  ambient  temperature 
of  25  °C  for  temperature  relaxation  to  the  environmental  temper¬ 
ature.  By  comparing  the  experimental  data  with  the  simulated  re¬ 
sults,  the  lumped  heat  transfer  coefficient  h  is  determined  to  be 
2.3  W  m”2  K_1,  which  represents  the  synthetic  action  of  convection 
and  radiation  heat  transfer. 

The  initial  conditions  of  the  model  are: 

C1  =  cl,0?  C2  =  C2,(n  T  =  To,  01, pe  =  01,term+  01,0s 
01, ne  —  01, term-  =  0,  02, ne  =  02, sp  =  02, pe  =  02,0 


3.  Numerical  simulation  and  experimental  validation 

3.1.  Numerical  simulation 

A  finite  element  based  commercial  software,  COMSOL  Multi- 
physics  4.3,  was  employed  to  solve  the  2D  electrochemical-thermal 
coupled  model  presented  in  this  paper.  A  given  battery  current  iapp 
was  applied  at  the  terminal  boundary.  The  outputs  of  the  model  are 
the  battery  working  voltage,  electronic  and  ionic  current  density 
distribution,  lithium  and  li-ion  concentration  distribution  and  tem¬ 
perature  field.  The  MUMPS  time  dependent  solver  was  chosen  with  a 
relative  tolerance  of  10-3  for  all  the  variables;  solutions  were  tested 
for  mesh  independence;  and  the  solution  time  for  this  model  was 
within  45  min  by  the  computing  platform  with  quad-core  processors 
(3.2  GHz)  and  a  total  of  8  GB  random  access  memory  (RAM). 

3.2.  Experimental  validation 

A  commercial  Type  38120  LiFeP04  power  battery  (</> 
38  mm  x  120  mm  cylindrical  battery,  capacity  10  Ah,  nominal 
voltage  3.2  V,  1.5  mol  L-1  LiPF6  in  EC/DMC  solvent)  was  experi¬ 
mentally  tested.  A  YKLE  test  system  was  used  to  monitor  the 
charge-discharge  current  and  working  voltage.  As  shown  in  Fig.  4, 


At  the  outer  can  surface  V,  the  lumped  heat  exchange  between 
battery  and  environment  is  described  in  Eq.  (30) 


Fig.  3.  Determination  of  lumped  heat  transfer  coefficient  h  by  comparing  experimental 
data  with  simulated  results. 
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Fig.  4.  Top  view  of  the  battery  ceil  with  the  placement  of  the  thermocouples. 


six  thermocouples  (from  negative  terminal  to  positive  terminal 
indexed  as  #1  through  #6)  were  evenly  placed  on  the  side  face  of  the 
battery  within  the  same  axial  section  to  monitor  the  surface  tem¬ 
perature  change.  According  to  the  manufacturer,  the  end-of- 
discharge  voltage  is  2.6  V  at  constant  current  discharge;  the  stop 
condition  for  constant  charge  process  is  that  being  charged  to  3.65  V 
and  then  with  constant  voltage  at  3.65  V  till  the  charge  current 
declines  to  0.2  A.  Another  stop  condition  for  both  charge  and 
discharge  is  when  the  total  charge/discharge  capacity  reaches  10  Ah. 

The  model  validation  is  conducted  by  comparing  simulation 
results  with  experimental  data  of  charging  and  discharging  be¬ 
haviors  at  different  conditions.  The  focus  will  be  on  thermal  vali¬ 
dation  and  electrochemical  validation. 


3.2.1.  Thermal  validation 

Fig.  5  shows  the  average  temperature  profile  of  the  battery 
surface  obtained  by  both  simulation  and  experiment  measurement. 
The  operating  condition  for  the  battery  in  Fig.  5  is  to  discharge  the 
battery  under  different  C-rates  (0.3C,  0.5C,  1C,  2C)  at  constant 
ambient  temperature  (25  °C).  The  average  measured  surface  tem¬ 
perature  is  an  average  of  temperatures  obtained  by  the  six  ther¬ 
mocouples.  As  shown  in  Fig.  5,  the  average  surface  temperature 
increases  by  around  1.4  K,  4  K,  8.1  K,  and  15.8  K  above  the  ambient 
temperature  for  discharge  rates  of  0.3C,  0.5C,  1C,  and  2C,  respec¬ 
tively.  It  also  can  be  seen  that  the  surface  temperature  have  the 
similar  growth  trends  for  the  low  discharge  rates,  while  the  tem¬ 
perature  increases  very  fast  at  a  high  discharge  rate,  in  particular 
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Fig.  6.  Comparison  of  simulated  results  of  axial  surface  temperatures  with  experi¬ 
mental  data  during  1C  discharge  under  natural  convection  condition. 


when  the  discharge  rate  is  above  1 C  in  this  study.  The  variation  of 
temperature  profile  may  be  due  to  the  internal  heat  generation, 
which  will  be  discussed  in  Section  4.  Fig.  6  shows  the  axial  tem¬ 
perature  profiles  of  the  battery  surface  obtained  by  both  simulation 
and  measurement.  The  operating  condition  for  the  battery  in  Fig.  6 
is  to  discharge  the  battery  under  1C  rate  at  constant  ambient 
temperature  (25  °C).  The  surface  temperature  was  measured  by  six 
thermocouples  evenly  attached  to  the  surface  of  the  battery  as 
shown  in  Fig.  4.  The  simulated  results  differ  slightly  from  the 
measured  data  towards  the  end  of  discharge.  The  simulated  tem¬ 
perature  appears  low  in  the  regions  near  both  terminals,  while  the 
experimental  data  shows  that  these  regions  have  higher  tempera¬ 
ture.  This  is  because  the  electrical  contact  resistance  between  the 
connectors  and  terminals  of  the  cell  would  cause  temperature  rise 
due  to  the  local  ohmic  heat  generation  [12].  The  contact  resistance, 
however,  is  not  included  in  the  model. 

3.2.2.  Electrochemical  validation 

Fig.  7  shows  the  battery  working  voltage  during  discharge 
process  under  different  C-rates  (0.3C,  1C,  2C)  at  constant  ambient 
temperature  (25  °C).  The  simulation  results  agree  with  the 


Discharge  capacity  (Ah) 

Fig.  5.  Comparison  of  simulated  results  of  average  surface  temperatures  with  exper¬ 
imental  data  during  galvanostatic  discharge  (0.3C,  0.5C,  1C,  2C)  under  natural  con¬ 
vection  condition. 


Discharge  capacity  (Ah) 


Fig.  7.  Comparison  of  simulated  results  of  working  voltage  with  experimental  data 
during  galvanostatic  discharge  (0.3C,  1C,  2C)  under  natural  convection  condition. 
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Fig.  8.  Schematic  diagram  of  current  distribution  in  a  local  cell  unit. 


experimental  data  generally  well  over  the  entire  discharge  process. 
Yet,  there  are  a  few  differences  in  general  discharge  voltage  plateau 
especially  towards  the  end  of  discharge  cycle.  It  is  important  to 
emphasize  that  these  differences  are  related  to  the  thermodynamic 
data  and  the  design  of  batteries.  For  example,  the  half-cell  open 
circuit  potentials  for  both  LixFeP04  and  LixC6  electrodes  used  in  this 
model  are  not  measured  values  for  the  experimental  battery  used 
here  but  instead  from  published  literature;  the  actual  capacity  for  a 
new  battery  is  always  lager  than  that  offered  by  the  manufacturer, 
which  has  an  impact  on  the  initial  battery  capacity  design.  If  these 
data  were  made  available  for  the  same  battery  used  here  for  model 
validation,  the  remaining  discrepancy  between  the  experimental 
data  and  simulation  results  would  have  been  further  reduced. 


4.  Results  and  discussion 

4.1.  Electrochemical  performance 

Fig.  8  shows  the  current  flow  in  the  local  cell  unit  of  the  battery. 
In  the  discharging  process,  the  current  in  the  x  direction  within  the 
electrodes  and  separator  is  so  small  that  the  current  flow  between 
the  two  current  collectors  is  assumed  to  be  perpendicular  to  the 
electrodes.  i0  represents  the  current  density  transferred  through 
the  separator  from  the  negative  electrode  to  the  positive  electrode, 
which  can  be  calculated  by  Eqs.  (7)  and  (8)  subject  to  the  associated 
boundary  conditions  (23)— (24).  We  know  that  the  current  density 
io  transferring  through  the  local  cell  unit  is  affected  strongly  by  the 
potential  of  both  electrodes  which  can  be  calculated  by  the  current 
continuity  equation. 

Fig.  9  shows  the  potential  distribution  on  the  positive  and 
negative  electrodes  during  the  discharge  with  1C  rate.  It  can  be  seen 
that  the  values  of  the  potential  of  Fig.  9(a)— (c)  are  around  3.23, 3.19, 
and  2.78  V  respectively.  By  halfway  discharge,  the  potential  drops 
slightly  complying  with  the  stable  discharge  voltage  platform  of 
LiFeP04  battery.  Flowever,  the  values  of  the  potential  and  the 
appearance  of  the  potential  distribution  on  the  negative  current 
collectors  during  the  discharge  has  hardly  changed,  thus  only  the 
potential  distribution  on  the  negative  current  collectors  at  the 
discharge  time  of  3  min  is  shown  in  Fig.  9(d).  This  is  due  to  the  fact 
that  the  potential  of  the  negative  terminal  is  specially  defined  as  0  V 
in  the  model.  Besides,  the  potentials  on  both  positive  and  negative 
current  collectors  have  the  similar  variation  trend  along  the  axial 
direction  of  the  battery.  As  shown  in  Fig.  8,  all  the  currents  inside  the 
battery  flow  from  the  negative  terminal  to  the  positive  terminal 
during  the  discharge,  which  would  suggest  that  the  values  of  the 
potential  decrease  along  the  current  direction  as  in  Fig.  9.  Therefore, 
in  a  cylindrical  battery,  the  potential  difference  between  positive 
and  negative  current  collectors  in  the  local  region  of  the  battery  axial 
cross-section  is  relative  uniform  as  compared  to  that  of  the  prismatic 
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Fig.  9.  Potential  distributions  on  the  positive  current  collector  (a,  b,  c)  and  the  negative  current  collector  (d)  during  discharge  with  1C  rate.  Discharge  time  is:  (a,  d)  3  min,  (b) 
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Fig.  10.  Distributions  of  current  density  transferred  through  the  separator  in  the  axial  cross-section  of  battery  during  discharge  at  1C  rate.  Discharge  time  is:  (a)  3  min,  (b)  30  min, 
and  (c)  57  min. 


battery  published  in  the  literature  1].  This  is  because  the  positive 
and  negative  terminals  are  placed  in  the  same  side  of  the  prismatic 
battery,  whereas  the  directions  of  the  current  on  the  positive  and 
negative  current  collectors  are  completely  opposite,  thus  leading  to 
the  non-uniform  distributions  of  local  potential  difference  as 
compared  to  that  of  the  cylindrical  battery  here. 

Inside  the  cell  unit,  the  ionic  current  transferring  through  the 
separator  is  equal  to  the  total  current  in  the  electrode  layers. 
Therefore  the  magnitude  of  the  ionic  current  transferring  through 
the  separator  reflects  the  reaction  rate.  Fig.  10  shows  the  non- 
uniform  distributions  of  current  density  transferred  through  the 
separator  in  the  axial  cross-section  of  battery  under  1C  rate 
discharge.  At  the  beginning  of  discharge  (Fig.  10(a)),  the  current 
density  decreases  from  the  negative  tabs  to  the  positive  ones  in  the 


Time  (s) 


Fig.  11.  Time  history  of  the  total  heat  generation  rate  in  the  axial  cross-section  of 
battery  regardless  of  the  current  collectors  during  discharge  at  1C  rate. 


axial  direction  y  coordinate.  Meanwhile,  in  the  radial  direction,  the 
current  density  also  decreases  in  x  coordinate.  Nevertheless,  after 
halfway  through  discharge  (Fig.  10(b)  and  (c)),  the  region  with 
higher  current  density  gradually  shifts  from  near  the  negative  tabs 
and  the  battery  axis  to  near  the  positive  tabs  and  the  outer  can  of 
the  battery.  This  is  because  the  region  with  the  higher  current 
density  at  the  initial  stage  becomes  more  utilized  by  the  halfway 
discharge  and  the  region  less  utilized  tends  to  have  a  higher  reac¬ 
tion  rate. 

4.2.  Heat  generation  and  temperature  distribution 

The  rise  in  temperature  within  a  cell  unit  originates  from  the 
electrochemical  reaction,  ohmic  heat  and  polarization  heat  as  the 


Fig.  12.  Ohmic  heat  distributions  on  current  collectors  during  discharge  at  1C  rate. 
Discharge  time:  3  min. 
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Fig.  13.  Distribution  of  temperature  at  the  end  of  discharge  at  (a)  1C,  (b)  5C  and  (c)  10C  rate  under  natural  convection  condition. 


charge-discharge  is  passed  through  the  battery.  Fig.  11  plots  the 
total  heat  generation  rates  under  1C  discharge  within  the  cell  units 
at  three  different  regions  of  Ai,  A2  and  A3  as  labeled  in  Fig.  1.  The 
heat  generation  rate  varies  very  little  from  each  other  among  the 
different  regions  during  the  entire  discharge  cycle.  The  region  Ai 
close  to  the  negative  tabs  seems  to  have  the  higher  heat  generation 
rate  than  do  the  region  A2  in  the  middle  and  the  region  A3  close  to 
the  positive  tabs  at  the  beginning  of  discharge,  but  this  trend  re¬ 
verses  after  the  halfway  through  discharge.  The  reason  is  that  the 


heat  generation  at  the  beginning  of  discharge  is  mainly  caused  by 
the  activation  reaction  rate  (working  current  density).  By 
comparing  Fig.  10  with  Fig.  11,  we  can  see  that  the  higher  reaction 
rate  the  local  cell  unit  has,  the  more  quickly  the  heat  is  generated  in 
the  local  region.  The  region  with  the  highest  working  current 
density  shifts  from  the  region  Ai  close  to  the  negative  tabs  to  the 
region  A3  close  to  the  positive  tabs  during  the  discharge.  This  trend 
agrees  with  that  of  the  heat  generation  rate.  Fig.  12  plots  the  ohmic 
heat  distributions  on  the  positive  and  negative  current  collectors 
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along  the  axial  direction.  It  can  be  seen  that  the  ohmic  heat 
generated  on  both  the  positive  and  the  negative  current  collectors 
decreases  from  one  end  near  the  tab  to  the  other  end,  which  is 
consistent  with  the  flow  of  electron  current  as  shown  in  Fig.  8.  In 
addition,  the  conductivity  of  the  positive  current  collector  made  by 
Al  foil,  is  lower  than  that  of  the  negative  current  collector  made  by 
Cu,  which  results  in  more  ohmic  heat  generated  on  the  positive 
current  collectors.  However,  it  is  also  of  interest  to  see  that  the 
summation  of  the  ohmic  heat  in  both  current  collectors  along  axial 
direction  is  relative  uniform  as  compared  to  that  in  the  lengthwise 
direction  of  the  prismatic  battery  in  literature  [2].  This  corresponds 
to  the  fact  that  the  positive  and  negative  tabs  are  placed  at  both 
ends  of  the  cylindrical  battery,  which  has  a  significant  effect  on  the 
current  density  and  potential  distributions,  thus  affecting  the 
temperature  distribution  of  the  battery. 

Fig.  13  shows  the  temperature  distribution  of  the  axial  cross- 
section  of  battery  during  discharge  at  1C,  5C  and  IOC  rate  under 
natural  convection  condition.  The  maximum  temperature  differ¬ 
ence,  between  the  hottest  and  the  coldest  spots,  at  1C,  5C  and  IOC 
rates  are  0.5  K,  3.1  K,  and  9.3  K,  respectively,  at  the  end  of  discharge; 
the  temperature  gradient  is  higher  in  the  axial  (x)  direction  than  that 
in  the  radial  (y )  direction.  The  reason  is  that  the  heat  dissipates  much 
faster  in  the  region  close  to  the  battery  terminal  sides  than  in  the 
battery  core  region,  since  the  region  near  both  terminal  sides  has 
relative  large  dissipated  surface  compared  to  the  region  close  to  the 
battery  core.  Fig.  13  also  shows  that  the  temperature  is  distributed 
symmetrically  around  the  cross-section  in  the  middle  of  the  axis  at 
all  three  discharge  rates.  This  phenomenon  can  be  further  explained 
by  the  results  shown  in  Figs.  11  and  12.  As  indicated  by  Fig.  11,  the 
overall  heat  generate  rates  vary  very  little  along  the  axial  direction. 
Also  the  total  ohmic  heat  generated  on  the  current  collectors  are 
relative  uniform  in  the  axial  direction  as  shown  in  Fig.  12.  The 
positioning  of  the  current  collectors  is  one  of  the  primary  causes  to 
have  such  a  symmetric  temperature  distribution  since  the  current 
collector  determines  the  potential  distribution,  which  further  af¬ 
fects  the  current  distribution  and  the  heat  generation. 

5.  Conclusion 

A  two-dimensional  electrochemical— thermal  model  has  been 
developed  for  a  commercial  Type  38120  LiFeP04  cylindrical  power 
battery  by  coupling  the  mass,  charge,  and  energy  conservations  as 
well  as  the  electrochemical  kinetics.  The  effect  of  the  current  col¬ 
lecting  tabs  on  the  battery  has  been  taken  into  consideration.  The 
modeling  results  are  further  validated  for  both  electrochemical 
performances  and  thermal  behavior  using  experimental  data.  The 
modeling  results  indicate  that  the  potentials  on  both  positive  and 
negative  current  collectors  have  the  similar  variation  trend  along 
the  axial  direction  of  the  battery  due  to  the  arrangement  of  the 
positive  and  negative  current  collecting  tabs  in  the  cylindrical 
battery.  The  distribution  of  potential  difference  in  the  axial  cross- 
section  of  the  cylindrical  battery  is  much  more  uniform  than  that 
in  the  lengthwise  direction  of  the  prismatic  battery. 

The  reaction  rate  and  the  total  heat  generation  rate  of  the 
electrodes  are  also  evaluated.  The  higher  reaction  rate  the  local  cell 
unit  has,  the  more  quickly  the  heat  is  generated  in  the  local  region. 
The  heat  generation  rate  among  the  different  regions  is  relative 
uniform  during  the  entire  discharge  cycle. 

Compared  to  a  prismatic  battery  in  the  lengthwise  direction,  the 
distribution  of  the  sum  of  the  ohmic  heat  generated  on  the  positive 
and  negative  current  collectors  is  relative  uniform  along  the  axial 
direction  of  the  cylindrical  battery. 

The  positioning  of  the  current  collectors  is  one  of  the  primary 
causes  to  have  such  a  symmetric  temperature  distribution  on  the 
cylindrical  battery  since  the  current  collector  determines  the 


potential  distribution,  which  further  affects  the  current  distribution 
and  the  heat  generation. 
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Nomenclature 

Ci  concentration  of  lithium  in  the  active  material  particles 

(mol  m-3) 

C2  electrolyte  concentration  (mol  m-3) 

Ci, max  maximum  concentration  of  lithium  in  the  active  material 
(mol  m-3) 

Ci, surf  surface  concentration  of  lithium  in  the  active  material 
(mol  m-3) 

Di  diffusion  coefficient  of  lithium  in  the  active  material 

(m2  s-1) 

D2  diffusion  coefficient  of  electrolyte  (m2  s-1) 

Eao  diffusion  coefficient  activation  energy  (kj  mol-1) 

EaR  reaction  activation  energy  (kj  mol-1) 

Eceii  working  voltage  of  the  battery  (V) 
f±  average  molar  activity  coefficient 

E  Faraday’s  constant  (C  mol-1) 

h  lumped  heat  transfer  coefficient  (W  m-2  I<-1) 

iapp  applied  current  density  of  the  battery  (A  m-2) 
jo  exchange  current  density  (A  m-2) 

jn  local  charge  transfer  current  density  (A  m-2) 

ko  reaction  rate  constant  (m2,5  mol-0,5  s-1) 

1<  thermal  conductivity  (W  (m  I<)-1) 

q  volumetric  heat  generation  (W  m-3) 

qrea  reaction  heat  generation  (W  m-3) 
qac t  polarization  heat  generation  (W  m-3) 
qohm  ohmic  heat  generation  (W  m-3) 

R  gas  constant,  8.314  (J  mol-2  I<-1) 

r  radius  distance  variable  of  electrode  particles  (pm) 

ro  radius  of  electrode  particles  (pm) 

Sa  specific  surface  area  (m-1) 

SOC  state  of  charge 

t  time  (s) 

t+  transferring  number  of  Li+ 

T  absolute  temperature  (I<) 

Tamb  ambient  temperature  (I<) 

Ueq  open  circuit  potential  of  the  electrode  (V) 

Uref  open  circuit  potential  under  the  reference  temperature 

(V) 

Greek  letters 

aa  anodic  transfer  coefficient 

ac  cathodic  transfer  coefficient 

e\  active  material  volume  fraction 

£2  electrolyte  volume  fraction 

0i  solid  phase  potential  (V) 

02  liquid  phase  potential  (V) 

y  Bruggeman  tortuosity  exponent 

g\  electronic  conductivity  of  solid  matrix  (S  m-1) 

(72  ionic  conductivity  of  electrolyte  (S  m-1) 

5 i  thickness  of  each  battery  component  (pm) 

rj  local  surface  overpotential  (V) 

Subscripts ,  superscripts  and  acronyms 
0  initial  or  equilibrated  value 

1  solid  phase 
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2  liquid  phase 

amb  ambient 

cc  current  collector 

eff  effective  value 

el  electrolyte 

ne  negative  electrode 

pe  positive  electrode 

ref  reference  value 

sp  separator 

surf  surface  of  active  material  particles 

term +/-  positive/negative  terminal 
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